#include "Modeling\CUDADataStruct.h"
//#include ".\Modeling\RenderBaseStruct.h"
#include ".\CUDA\utility\inc\cutil.h"
#include ".\CUDA\utility\inc\cudpp\cudpp.h"
#include <stdio.h>

typedef unsigned int uint;

#define UNIFORM_KNOT 1
#define NORMALFUNCTION 0

#define NTHREADS		128
#define PIXELTHREADS	64
//const int nThreads = 128;


#define DEBUG_BUFFERS 0
extern void dumpBuffer(uint *d_buffer, int nelements);

const uint pixelBits = 10;
uint patchBits, pixelMask, patchMask;


uint3 deg, degShift, seg, segShift;
float3 boxSize, box0;

uint nPatch;

#ifndef UNIFORM_KNOT
__constant__ float* dBoxX = 0;
__constant__ float* dBoxY = 0;
__constant__ float* dBoxZ = 0;
#endif

__constant__ int* dPatchIndex = 0;
__constant__ float* dFuncData = 0;


//the variable below are defined in the way described in "Progress 06.22"
uint maxSeg;
uint maxSilPoint;
//uint maxPolarIndexBlock;
//uint nPatchEveryBlock;


size_t* dMaxSilPoint  = NULL;
__device__ uint* dSilOccupied = 0;
__device__ uint* dSilOccupiedCompact = 0;
__device__ uint* dSilOccupiedValid = 0;
__device__ uint* dSilOccupiedSort = 0;
__device__ uint* dSilOccupiedSortPad = 0;
//__device__ uint* dSilOccupiedSortIndex = 0;

uint* polarOccupiedSort;
uint* polarOccupiedPad;

CUDPPHandle plan;

__constant__ MULTI_SIL_RENDER* dMultiSilRender;
int renderWidth, renderHeight;
//float4 viewProj0, viewProj1, viewProj2, viewProj3;
float4 invViewProj0, invViewProj1, invViewProj2, invViewProj3;
//float3 eye;
//float errorXTol;
//LIGHT light[3];
//MATERIAL material;
//int nLight;

__device__ void normalize(float* a)
{
	float sum;
	sum=sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]);
	if(sum > 1.0e-10f)
    {
		a[0]/=sum; a[1]/=sum; a[2]/=sum;
	}
}

__device__ void addambient(MATERIAL *mat, float rgb[])
{
	rgb[0]+=mat->ka*mat->d[0]*mat->a[0];
	rgb[1]+=mat->ka*mat->d[1]*mat->a[1];
	rgb[2]+=mat->ka*mat->d[2]*mat->a[2];
}

__device__ void addPhong(MATERIAL *mat, LIGHT* light, float norm[], float ray[], float rgb[])
{
	float H[3] ={(-light->dir[0]-ray[0])/2.0f, (-light->dir[1]-ray[1])/2.0f, (-light->dir[2]-ray[2])/2.0f};
	normalize(H);
	float phong = norm[0]*H[0]+norm[1]*H[1]+norm[2]*H[2];
	phong *= phong;
    phong *= phong;
    phong *= mat->ks;
	rgb[0]+=phong*light->rgb[0];
	rgb[1]+=phong*light->rgb[1];
	rgb[2]+=phong*light->rgb[2];
}

__device__ void addDiffuse(MATERIAL *mat, LIGHT* light, float diff, float norm[], float rgb[])
{
	rgb[0]+=mat->kd*diff*light->rgb[0];
	rgb[1]+=mat->kd*diff*light->rgb[1];
	rgb[2]+=mat->kd*diff*light->rgb[2];
}

__device__ void GaussEliminationGPU(float A[][3], float b[])
{
	int i,j,k,ik;
	float mik, temp;
	for(k=0; k<3; k++)
	{
		mik=-1.0f;
		for(i=k; i<3; i++)
		{
			if(fabs(A[i][k])>mik)
			{
				mik=fabs(A[i][k]);
				ik=i;
			}
		}
		if (mik<1.0e-10f)
		{
			b[0]=b[1]=b[2]=0.0f;
			return;
		}
		for(j=k; j<3; j++)
		{
			//swap(A[ik][j], A[k][j]);
			temp = A[ik][j];
			A[ik][j]=A[k][j];
			A[k][j]=temp;
		}
		//swap(b[ik], b[k]);
		temp = b[ik];
		b[ik] = b[k];
		b[k] = temp;
		b[k] /= A[k][k];
		for(i=3-1; i>=k; i--)
			A[k][i]/=A[k][k];
		for(i=k+1; i<3; i++)
		{
			b[i] -= A[i][k]*b[k];
			for (j=3-1; j>k; j--)
				A[i][j] -= A[i][k]*A[k][j];
		}
	}
	for(i=3-1; i>=0; i--)
	{
		for (j=i+1; j<3; j++)
			b[i] -= A[i][j]*b[j];
	}
}

__device__ void GetBernstein(float Bern[], int deg, float u)
{
	float cu = 1.0f-u;
	switch (deg)
	{
	case 2:
		Bern[0]=cu*cu; Bern[1]=2*u*cu; Bern[2]=u*u; break;
	case 1:
		Bern[0]=cu; Bern[1]=u; break;
	case 3:
		Bern[0]=cu*cu*cu; Bern[1]=3*u*cu*cu; Bern[2]=3*u*u*cu; Bern[3]=u*u*u; break;
	}
}

__device__ void EvalBezierBlossom(float u, float v, float w, int nx, int ny, int nz, float* func, float* re)
{
	float nn[4];
	float f1[4][4];
	float f2[4];
	GetBernstein(nn, nz, w);
	for (int i=0; i<=nx; i++)
	for (int j=0; j<=ny; j++)
	{
		f1[i][j] = 0.0f;
		for (int k=0; k<=nz; k++)
			f1[i][j] += *func++ * nn[k];
	}

	GetBernstein(nn, ny, v);
	for (int i=0; i<=nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<=ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}

	GetBernstein(nn, nx, u);
	(*re) = 0.0f;
	for (int i=0; i<= nx; i++)
	{
		(*re) += f2[i]*nn[i];
	}
}

__device__ void EvalBezierBlossomNormal(float u, float v, float w, int nx, int ny, int nz, uint3 degShift, float* func, float* re)
{
	float nn[4];
	float f1[4][4];
	float f2[4];

	//X direction:
	GetBernstein(nn, nz, w);
	for (int i=0; i<nx; i++)
	for (int j=0; j<=ny; j++)
	{
		f1[i][j] = 0.0f;
		for (int k=0; k<=nz; k++)
			f1[i][j] += *func++ * nn[k];
	}
	func += degShift.x;
	GetBernstein(nn, ny, v);
	for (int i=0; i<nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<=ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}
	GetBernstein(nn, nx-1, u);
	re[0] = 0.0f;
	for (int i=0; i<nx; i++)
	{
		re[0] += f2[i]*nn[i];
	}

	//X direction:
	GetBernstein(nn, nz, w);
	for (int i=0; i<=nx; i++)
	{
		for (int j=0; j<ny; j++)
		{
			f1[i][j] = 0.0f;
			for (int k=0; k<=nz; k++)
				f1[i][j] += *func++ * nn[k];
		}
		func += degShift.y;
	}
	GetBernstein(nn, ny-1, v);
	for (int i=0; i<=nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}
	GetBernstein(nn, nx, u);
	re[1] = 0.0f;
	for (int i=0; i<= nx; i++)
	{
		re[1] += f2[i]*nn[i];
	}

	//Z direction:
	GetBernstein(nn, nz-1, w);
	for (int i=0; i<=nx; i++)
	{
		for (int j=0; j<=ny; j++)
		{
			f1[i][j] = 0.0f;
			for (int k=0; k<nz; k++)
				f1[i][j] += *func++ * nn[k];
			func++;
		}
	}
	GetBernstein(nn, ny, v);
	for (int i=0; i<=nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<=ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}
	GetBernstein(nn, nx, u);
	re[2] = 0.0f;
	for (int i=0; i<= nx; i++)
	{
		re[2] += f2[i]*nn[i];
	}
}

extern "C"
void InitMulitSilSurfaceCUDA(float* pFuncData, int* Index, int* indexInverse, int NumPatch, MULTI_SIL_SURFACE* multiSilSurface, int surfWidth, int surfHeight)
{
	//Attention: not done with indexInverse;
	patchBits = 8*sizeof(uint)-2*pixelBits;
	uint maxNumPatch = 1<<patchBits;
	if (NumPatch>maxNumPatch)
	{
		printf("Too many patchs!\n"); 
		return;
	}

	pixelMask = (1<<pixelBits)-1;
	patchMask = (1<<patchBits)-1;


	nPatch = NumPatch;
	deg.x=multiSilSurface->degX; deg.y=multiSilSurface->degY; deg.z=multiSilSurface->degZ;
	seg.x=multiSilSurface->segX; seg.y=multiSilSurface->segY; seg.z=multiSilSurface->segZ;
	degShift.x = (deg.y+1)*(deg.z+1); degShift.y = deg.z+1; degShift.z = 1;
	segShift.x = seg.y*seg.z; segShift.y = seg.z; segShift.z = 1;
	maxSeg = max(seg.x, seg.y); maxSeg = max(maxSeg, seg.z);

	int funcSize = NumPatch*(deg.x+1)*(deg.y+1)*(deg.z+1)*sizeof(float);	
	CUDA_SAFE_CALL( cudaMalloc((void**)&dFuncData, funcSize) );
	CUDA_SAFE_CALL( cudaMemcpy(dFuncData, pFuncData, funcSize, cudaMemcpyHostToDevice) );
	//cudaMemcpyToSymbol(dFuncData, pFuncData, size);	

	CUDA_SAFE_CALL( cudaMalloc((void**)&dPatchIndex, seg.x*seg.y*seg.z*sizeof(int)) );
	CUDA_SAFE_CALL( cudaMemcpy(dPatchIndex, Index, seg.x*seg.y*seg.z*sizeof(int), cudaMemcpyHostToDevice) );
	//cudaMemcpyToSymbol(dPatchIndex, Index, segX*segY*segZ*sizeof(uint));


	CUDA_SAFE_CALL( cudaMalloc((void**)&dMultiSilRender, sizeof(MULTI_SIL_RENDER)) );
	renderWidth = surfWidth;
	renderHeight = surfHeight;
 	int silSize = renderWidth*renderHeight*sizeof(uint);
	CUDA_SAFE_CALL( cudaMalloc((void**)&dSilOccupied, silSize) );		
	CUDA_SAFE_CALL( cudaMalloc((void**)&dSilOccupiedValid, silSize) );
	//maxSilPoint = renderHeight*renderWidth;
	CUDA_SAFE_CALL( cudaMalloc((void**)&dSilOccupiedCompact, silSize) );
	CUDA_SAFE_CALL( cudaMalloc((void**)&dSilOccupiedSort, silSize) );
    CUDA_SAFE_CALL( cudaMalloc( (void**) &dMaxSilPoint, sizeof(size_t)));

	//uint numPolarPad = silSize*PIXELTHREADS*3/5;
	uint numPolarPad = silSize+sizeof(uint);
	polarOccupiedSort = (uint*)malloc(numPolarPad);
	polarOccupiedPad = (uint*)malloc(numPolarPad*4/3);
	printf("Attention: We got memory leak here!\n");
	CUDA_SAFE_CALL( cudaMalloc((void**)&dSilOccupiedSortPad, numPolarPad*4/3) );


#if	UNIFORM_KNOT	
	box0.x = multiSilSurface->boxX[0];
	box0.y = multiSilSurface->boxY[0];
	box0.z = multiSilSurface->boxZ[0];
	boxSize.x = multiSilSurface->boxX[1]-multiSilSurface->boxX[0];
	boxSize.y = multiSilSurface->boxY[1]-multiSilSurface->boxY[0];
	boxSize.z = multiSilSurface->boxZ[1]-multiSilSurface->boxZ[0];
#else	
	CUDA_SAFE_CALL( cudaMalloc((void**)&dBoxX, (seg.x+1)*sizeof(float)) );
	CUDA_SAFE_CALL( cudaMalloc((void**)&dBoxY, (seg.y+1)*sizeof(float)) );
	CUDA_SAFE_CALL( cudaMalloc((void**)&dBoxZ, (seg.z+1)*sizeof(float)) );
	CUDA_SAFE_CALL( cudaMemcpy(dBoxX, multiSilSurface->boxX, (seg.x+1)*sizeof(float), cudaMemcpyHostToDevice) );	
	CUDA_SAFE_CALL( cudaMemcpy(dBoxY, multiSilSurface->boxY, (seg.y+1)*sizeof(float), cudaMemcpyHostToDevice) );
	CUDA_SAFE_CALL( cudaMemcpy(dBoxZ, multiSilSurface->boxZ, (seg.z+1)*sizeof(float), cudaMemcpyHostToDevice) );
	//cudaMemcpyToSymbol(dBoxX, multiSilSurface->boxX, (seg.x+1)*sizeof(float));
	//cudaMemcpyToSymbol(dBoxY, multiSilSurface->boxY, (seg.y+1)*sizeof(float));
	//cudaMemcpyToSymbol(dBoxZ, multiSilSurface->boxZ, (seg.z+1)*sizeof(float));
#endif
}

__global__ void CheckPointUniform(unsigned char* surface, size_t pitch, int width, int height, float3 box0, float3 boxSize, uint3 seg, uint3 segShift, uint pixelBits, uint patchBits, int* dPatchIndex, uint* dSilOccupied, uint* dSilOccupiedValid, float4 invViewProj0, float4 invViewProj1, float4 invViewProj2, float4 invViewProj3)
{	
	uint x = blockIdx.x*blockDim.x + threadIdx.x;
	uint y = blockIdx.y*blockDim.y + threadIdx.y;
	if(x>=width || y>=height)		return;
	uint pixelI = y*width+x;
    float* pixel = (float*)(surface + y*pitch) + 4*x;
	//if( pixel[0] == 0.0f || pixel[0] == 1.0f)		return;
	if( pixel[0]<0.1f)								return;

	float posx = 2.0f*(x+0.5f)/width-1.0f;
	float posy = 1.0f-2.0f*(y+0.5f)/height;
	float posz = pixel[2];
	float worldW =  posx*invViewProj0.w + posy*invViewProj1.w + posz*invViewProj2.w + invViewProj3.w;
	float worldX = (posx*invViewProj0.x + posy*invViewProj1.x + posz*invViewProj2.x + invViewProj3.x) / worldW;
	float worldY = (posx*invViewProj0.y + posy*invViewProj1.y + posz*invViewProj2.y + invViewProj3.y) / worldW;
	float worldZ = (posx*invViewProj0.z + posy*invViewProj1.z + posz*invViewProj2.z + invViewProj3.z) / worldW;

	worldX = (worldX-box0.x)/boxSize.x;
	worldY = (worldY-box0.y)/boxSize.y;
	worldZ = (worldZ-box0.z)/boxSize.z;

	if(worldX<0 || worldX>=seg.x || worldY<0 || worldY>=seg.y || worldZ<0 || worldZ>=seg.z)
		return;
	else
	{
		//int patchI = dPatchIndex[(int)floorf(worldX)*segShift.x+(int)floorf(worldY)*segShift.y+(int)floorf(worldZ)*segShift.z];
		int patchI = (int)floorf(worldX)*segShift.x+(int)floorf(worldY)*segShift.y+(int)floorf(worldZ)*segShift.z;
		//if(patchI>0)
		{
			dSilOccupied[pixelI] = (((x<<pixelBits)|y)<<patchBits)|patchI;
			dSilOccupiedValid[pixelI] = 1;
		}
	}	
}

__global__ void UpdatePointKernel(unsigned char* surfaceIn, unsigned char* surfaceOut, size_t pitch, 
									 //int width, int height, float3 eye, float errorXTol, 
									 //float4 viewProj0, float4 viewProj1, float4 viewProj2, float4 viewProj3,
									 //float4 invViewProj0, float4 invViewProj1, float4 invViewProj2, float4 invViewProj3,
									 MULTI_SIL_RENDER* dMultiSilRender, //const uint allReadRenderNum, const uint otherReadRenderNum,
									 uint pixelBits, uint patchBits, 
									 //uint* SilOccupiedSort, uint*SilOccupiedSortIndex,
									 uint* SilOccupiedSortPad, float* FuncData, int* PatchIndex,
									 uint3 deg, uint3 degShift, uint3 seg, uint3 segShift, 
									 float3 boxSize, float3 box0, int iter)
{	
	//uint3 degShift;
	//degShift.z = 1;
	//degShift.y = deg.z+1;
	//degShift.x = (deg.y+1)*degShift.y;
	uint ind = threadIdx.x;
	__shared__ uint polarInfo[PIXELTHREADS];
	polarInfo[ind] = SilOccupiedSortPad[ind+blockIdx.x*blockDim.x];
	__syncthreads();

	uint polarIndex = (polarInfo[0]<<(pixelBits*2))>>(pixelBits*2);
	uint patchIndexX = polarIndex/segShift.x;
	uint patchIndexY = polarIndex%segShift.x/segShift.y;
	uint patchIndexZ = polarIndex%segShift.y;
	float patchX = patchIndexX*boxSize.x + box0.x;
	float patchY = patchIndexY*boxSize.y + box0.y;
	float patchZ = patchIndexZ*boxSize.z + box0.z;

	polarIndex = PatchIndex[polarIndex];

	extern __shared__ float funcNormal[];
	uint funcSize = (deg.x+1)*degShift.x;
	float* func = funcNormal;
#ifdef NORMALFUNCTION
	float* funcJacobi = &func[funcSize];
#else
#endif

	if(ind<funcSize)
		func[ind] = FuncData[polarIndex*funcSize+ind];
	__syncthreads();

#ifdef NORMALFUNCTION
	if((ind/degShift.x)<deg.x)
	{
		funcJacobi[ind]				= deg.x*(func[ind +degShift.x] - func[ind] );
	}
	if((ind%degShift.x/degShift.y)<deg.y)
	{
		funcJacobi[funcSize +ind]	= deg.y*(func[ind +degShift.y] - func[ind] );
	}
	if((ind%degShift.y)<deg.z)
	{
		funcJacobi[funcSize*2 +ind] = deg.z*(func[ind +degShift.z] - func[ind] );
	}
	__syncthreads();
#else
	float delta = 0.01f;
#endif

	if(polarInfo[ind]==0)
		return;


	float xyz[3] = {0.0f, 0.0f, 0.0f};
	float norm[3] = {0.0f, 0.0f, 0.0f};
	float fValue[3] = {0.0f, 0.0f, 0.0f};
	float horPlane[4] = {0.0f, 0.0f, 0.0f, 0.0f};
	float verPlane[4] = {0.0f, 0.0f, 0.0f, 0.0f};
	float jacobi[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};

	uint pixelX = (polarInfo[ind]>>(patchBits+pixelBits));
	uint pixelY = (polarInfo[ind]<<pixelBits)>>(patchBits+pixelBits);

	float screenX = 2.0f*(pixelX+0.5f)/dMultiSilRender->renderWidth-1.0f;
	float screenY = 1.0f-2.0f*(pixelY+0.5f)/dMultiSilRender->renderHeight;
	float screenZ = *((float*)(surfaceIn+pixelY*pitch) + pixelX*4 + 2);

	horPlane[0] = (screenY*dMultiSilRender->viewProj[0][3]-dMultiSilRender->viewProj[0][1]);
	horPlane[1] = (screenY*dMultiSilRender->viewProj[1][3]-dMultiSilRender->viewProj[1][1]);
	horPlane[2] = (screenY*dMultiSilRender->viewProj[2][3]-dMultiSilRender->viewProj[2][1]);
	horPlane[3] = (screenY*dMultiSilRender->viewProj[3][3]-dMultiSilRender->viewProj[3][1]) + patchX*horPlane[0] + patchY*horPlane[1] + patchZ*horPlane[2];
	horPlane[0] *= boxSize.x;
	horPlane[1] *= boxSize.y;
	horPlane[2] *= boxSize.z;
	verPlane[0] = (screenX*dMultiSilRender->viewProj[0][3]-dMultiSilRender->viewProj[0][0]);
	verPlane[1] = (screenX*dMultiSilRender->viewProj[1][3]-dMultiSilRender->viewProj[1][0]);
	verPlane[2] = (screenX*dMultiSilRender->viewProj[2][3]-dMultiSilRender->viewProj[2][0]);
	verPlane[3] = (screenX*dMultiSilRender->viewProj[3][3]-dMultiSilRender->viewProj[3][0]) + patchX*verPlane[0] + patchY*verPlane[1] + patchZ*verPlane[2];
	verPlane[0] *= boxSize.x;
	verPlane[1] *= boxSize.y;
	verPlane[2] *= boxSize.z;

	float worldW = screenX*dMultiSilRender->invViewProj[0][3] + screenY*dMultiSilRender->invViewProj[1][3] + screenZ*dMultiSilRender->invViewProj[2][3] + dMultiSilRender->invViewProj[3][3];
	xyz[0] = (screenX*dMultiSilRender->invViewProj[0][0] + screenY*dMultiSilRender->invViewProj[1][0] + screenZ*dMultiSilRender->invViewProj[2][0] + dMultiSilRender->invViewProj[3][0]) / worldW;
	xyz[1] = (screenX*dMultiSilRender->invViewProj[0][1] + screenY*dMultiSilRender->invViewProj[1][1] + screenZ*dMultiSilRender->invViewProj[2][1] + dMultiSilRender->invViewProj[3][1]) / worldW;
	xyz[2] = (screenX*dMultiSilRender->invViewProj[0][2] + screenY*dMultiSilRender->invViewProj[1][2] + screenZ*dMultiSilRender->invViewProj[2][2] + dMultiSilRender->invViewProj[3][2]) / worldW;

	xyz[0] = (xyz[0]-patchX)/boxSize.x;
	xyz[1] = (xyz[1]-patchY)/boxSize.y;
	xyz[2] = (xyz[2]-patchZ)/boxSize.z;

	float* pixel = (float*)(surfaceOut + pixelY*pitch) + 4*pixelX;

	float errF = 0.0f;
	float deltaPoint = 0.0f;
	for(int its=0; its<5; its++)
	{
		xyz[0] -= fValue[0];
		xyz[1] -= fValue[1];
		xyz[2] -= fValue[2];

		EvalBezierBlossom(xyz[0], xyz[1], xyz[2], deg.x, deg.y, deg.z, func, &fValue[0]);
		fValue[1] = horPlane[0]*xyz[0] + horPlane[1]*xyz[1] + horPlane[2]*xyz[2] + horPlane[3];
		fValue[2] = verPlane[0]*xyz[0] + verPlane[1]*xyz[1] + verPlane[2]*xyz[2] + verPlane[3];

#ifdef NORMALFUNCTION
		EvalBezierBlossomNormal(xyz[0], xyz[1], xyz[2], deg.x, deg.y, deg.z, degShift, funcJacobi, &jacobi[0][0]);
#else
		EvalBezierBlossom(xyz[0]+delta, xyz[1], xyz[2], deg.x, deg.y, deg.z, func, &jacobi[0][0]);
		EvalBezierBlossom(xyz[0], xyz[1]+delta, xyz[2], deg.x, deg.y, deg.z, func, &jacobi[0][1]);
		EvalBezierBlossom(xyz[0], xyz[1], xyz[2]+delta, deg.x, deg.y, deg.z, func, &jacobi[0][2]);
		jacobi[0][0] = (jacobi[0][0]-fValue[0])/delta;
		jacobi[0][1] = (jacobi[0][1]-fValue[0])/delta;
		jacobi[0][2] = (jacobi[0][2]-fValue[0])/delta;
#endif

		norm[0] = jacobi[0][0];
		norm[1] = jacobi[0][1];
		norm[2] = jacobi[0][2];
		jacobi[1][0] = horPlane[0];
		jacobi[1][1] = horPlane[1];
		jacobi[1][2] = horPlane[2];
		jacobi[2][0] = verPlane[0];
		jacobi[2][1] = verPlane[1];
		jacobi[2][2] = verPlane[2];

		errF = fValue[0]*fValue[0]+fValue[1]*fValue[1]+fValue[2]*fValue[2];
		if (errF <= dMultiSilRender->errorXTol*dMultiSilRender->errorXTol)		break;

		GaussEliminationGPU(jacobi, fValue);
		deltaPoint = fabsf(fValue[0])+fabsf(fValue[1])+fabsf(fValue[2]);
		if (deltaPoint <= dMultiSilRender->errorXTol)			break;	
	}


	MATERIAL material;
	material.ka = dMultiSilRender->material.ka;
	material.kd = dMultiSilRender->material.kd;
	material.ks = dMultiSilRender->material.ks;
	material.a[0] = dMultiSilRender->material.a[0];
	material.a[1] = dMultiSilRender->material.a[1];
	material.a[2] = dMultiSilRender->material.a[2];
	material.s[0] = dMultiSilRender->material.s[0];
	material.s[1] = dMultiSilRender->material.s[1];
	material.s[2] = dMultiSilRender->material.s[2];

	switch ((patchIndexX+patchIndexY+patchIndexZ)%4) 
	{
		case  0:
			material.d[0]=1.0; material.d[1]=material.d[2]=0.0;
			break;
		case  1:
			material.d[1]=1.0; material.d[0]=material.d[2]=0.0;
			break;
		case  3:
			material.d[2]=1.0; material.d[1]=material.d[0]=0.0;
			break;
		case  2:
			material.d[1]=0.0; material.d[2]=material.d[0]=1.0;
			break;
	}

    float rgb[3]={0.0f, 0.0f, 0.0f};

	norm[0] /= boxSize.x;
	norm[1] /= boxSize.y;
	norm[2] /= boxSize.z;
	normalize(norm);
	xyz[0] = xyz[0]*boxSize.x+patchX;
	xyz[1] = xyz[1]*boxSize.y+patchY;
	xyz[2] = xyz[2]*boxSize.z+patchZ;
    

	float ray[3]={xyz[0]-dMultiSilRender->eye[0], xyz[1]-dMultiSilRender->eye[1], xyz[2]-dMultiSilRender->eye[2]};
	normalize(ray);

    float diff;
	for(int i=0; i<dMultiSilRender->nLight; i++)  
	{
        diff = -(norm[0]*dMultiSilRender->light[i].dir[0]+norm[1]*dMultiSilRender->light[i].dir[1]+norm[2]*dMultiSilRender->light[i].dir[2]);
        if (diff>0)
        {
            //addPhong(&dMultiSilRender->material, &dMultiSilRender->light[i], norm, ray, rgb);
            //addDiffuse(&dMultiSilRender->material, &dMultiSilRender->light[i], diff, norm, rgb);
			addDiffuse(&material, &dMultiSilRender->light[i], diff, norm, rgb);
        }
	}
	//addambient(&dMultiSilRender->material, rgb);
	addambient(&material, rgb);
    pixel[0] = rgb[0];
    pixel[1] = rgb[1];
    pixel[2] = rgb[2];
}


extern "C"
void MultiSilRender(MULTI_SIL_RENDER* multiSilRender,	void* surfaceIn, void* surfaceOut, size_t pitch, int iter)
{
	double timeValue;
	unsigned int hTimer;
	CUT_SAFE_CALL(cutCreateTimer(&hTimer));
	printf("MultiSilRender ...\n");
	CUDA_SAFE_CALL(cudaThreadSynchronize());
	CUT_SAFE_CALL(cutResetTimer(hTimer));
	CUT_SAFE_CALL(cutStartTimer(hTimer));	

	CUDA_SAFE_CALL( cudaMemcpy(dMultiSilRender, multiSilRender, sizeof(MULTI_SIL_RENDER), cudaMemcpyHostToDevice) );
	//viewProj0.x = multiSilRender->viewProj[0][0]; viewProj0.y = multiSilRender->viewProj[0][1]; viewProj0.z = multiSilRender->viewProj[0][2]; viewProj0.w = multiSilRender->viewProj[0][3]; 
	//viewProj1.x = multiSilRender->viewProj[1][0]; viewProj1.y = multiSilRender->viewProj[1][1]; viewProj1.z = multiSilRender->viewProj[1][2]; viewProj1.w = multiSilRender->viewProj[1][3]; 
	//viewProj2.x = multiSilRender->viewProj[2][0]; viewProj2.y = multiSilRender->viewProj[2][1]; viewProj2.z = multiSilRender->viewProj[2][2]; viewProj2.w = multiSilRender->viewProj[2][3]; 
	//viewProj3.x = multiSilRender->viewProj[3][0]; viewProj3.y = multiSilRender->viewProj[3][1]; viewProj3.z = multiSilRender->viewProj[3][2]; viewProj3.w = multiSilRender->viewProj[3][3]; 
	invViewProj0.x = multiSilRender->invViewProj[0][0]; invViewProj0.y = multiSilRender->invViewProj[0][1]; invViewProj0.z = multiSilRender->invViewProj[0][2]; invViewProj0.w = multiSilRender->invViewProj[0][3]; 
	invViewProj1.x = multiSilRender->invViewProj[1][0]; invViewProj1.y = multiSilRender->invViewProj[1][1]; invViewProj1.z = multiSilRender->invViewProj[1][2]; invViewProj1.w = multiSilRender->invViewProj[1][3]; 
	invViewProj2.x = multiSilRender->invViewProj[2][0]; invViewProj2.y = multiSilRender->invViewProj[2][1]; invViewProj2.z = multiSilRender->invViewProj[2][2]; invViewProj2.w = multiSilRender->invViewProj[2][3]; 
	invViewProj3.x = multiSilRender->invViewProj[3][0]; invViewProj3.y = multiSilRender->invViewProj[3][1]; invViewProj3.z = multiSilRender->invViewProj[3][2]; invViewProj3.w = multiSilRender->invViewProj[3][3]; 
	//eye.x = multiSilRender->eye[0];
	//eye.y = multiSilRender->eye[1];
	//eye.z = multiSilRender->eye[2];
	//errorXTol = multiSilRender->errorXTol;
	//nLight = multiSilRender->nLight;
	//memset(&light, 0, sizeof(LIGHT)*3);
	//for(int i=0; i<nLight; i++)
	//{
	//	light[i].dir[0] = multiSilRender->light[i].dir[0];
	//	light[i].dir[1] = multiSilRender->light[i].dir[1];
	//	light[i].dir[2] = multiSilRender->light[i].dir[2];
	//	light[i].rgb[0] = multiSilRender->light[i].rgb[0];
	//	light[i].rgb[1] = multiSilRender->light[i].rgb[1];
	//	light[i].rgb[2] = multiSilRender->light[i].rgb[2];
	//}



 	int silSize = renderWidth*renderHeight*sizeof(uint);
	CUDA_SAFE_CALL( cudaMemset(dSilOccupied, 0, silSize) );
	CUDA_SAFE_CALL( cudaMemset(dSilOccupiedValid, 0, silSize) );

	//CollectSilPoint(surfaceIn, pitch, renderWidth, renderHeight, multiSilRender->invViewProj);	
	dim3 Db = dim3(16,16);
	dim3 Dg = dim3( (renderWidth+Db.x-1)/Db.x, (renderHeight+Db.y-1)/Db.y );

#if DEBUG_BUFFERS
	printf("dPatchIndex:\n");
	dumpBuffer((uint*)dPatchIndex, seg.x*seg.y*seg.z);
#endif	


#if	UNIFORM_KNOT	
	CheckPointUniform<<<Dg,Db>>>((unsigned char*)surfaceIn, pitch, renderWidth, renderHeight, 
								box0, boxSize, seg, segShift, pixelBits, patchBits,
								dPatchIndex, dSilOccupied, dSilOccupiedValid, 
								invViewProj0, invViewProj1, invViewProj2, invViewProj3);
#else
#endif

#if DEBUG_BUFFERS
	printf("dSilOccupied:\n");
	dumpBuffer(dSilOccupiedValid, renderWidth*renderHeight);
#endif

	//maxSilPoint = 3*renderWidth*maxSeg;
	//maxSilPoint = renderHeight*renderWidth;
	CUDA_SAFE_CALL( cudaMemset(dSilOccupiedCompact, 0, silSize) );

	CUDPPConfiguration config;
	CUDPPResult result = CUDPP_SUCCESS;	

    config.algorithm = CUDPP_COMPACT;
    config.datatype = CUDPP_UINT;
	config.options = CUDPP_OPTION_FORWARD;
	result = cudppPlan(&plan, config, renderWidth*renderHeight, 1, 0);
	if (result != CUDPP_SUCCESS)
	{
		printf("Error creating plan for compact\n");
		return;
	}
	result = cudppCompact(plan, dSilOccupiedCompact, dMaxSilPoint, dSilOccupied, dSilOccupiedValid, renderWidth*renderHeight);
	if (result != CUDPP_SUCCESS)
	{
		printf("Error compacting !\n");
		return;
	}	
	result = cudppDestroyPlan(plan);    
    if (result != CUDPP_SUCCESS)
    {
        printf("Error destroying CUDPPPlan!\n");
		return;
    }
	
	
	CUDA_SAFE_CALL( cudaMemcpy(&maxSilPoint, dMaxSilPoint, sizeof(size_t), cudaMemcpyDeviceToHost) );
#if DEBUG_BUFFERS
	printf("dSilOccupiedCompact:\n");
	dumpBuffer(dSilOccupiedCompact, maxSilPoint);
#endif



	CUDA_SAFE_CALL( cudaMemset(dSilOccupiedSort, 0, maxSilPoint*sizeof(uint)) );

	config.algorithm = CUDPP_SORT_RADIX_GLOBAL;
	config.datatype = CUDPP_UINT;
	config.op = CUDPP_ADD;
	config.options = (CUDPPOption)0;
	result = cudppPlan(&plan, config, maxSilPoint, 1, 0);
	if (result != CUDPP_SUCCESS)
	{
		printf("Error creating plan for SortBit\n");
		return;
	}
	result = cudppSortBit(plan, dSilOccupiedSort, dSilOccupiedCompact, maxSilPoint, patchBits);
	if (result != CUDPP_SUCCESS)
	{
		printf("Error SortBit\n");
		return;
	}
	result = cudppDestroyPlan(plan);    
    if (result != CUDPP_SUCCESS)
    {
        printf("Error destroying CUDPPPlan!\n");
    }


#if DEBUG_BUFFERS
	printf("dSilOccupiedSort:\n");
	dumpBuffer(dSilOccupiedSort, maxSilPoint);
#endif
	
	//GetSilSortIndex();

	CUDA_SAFE_CALL( cudaMemcpy(polarOccupiedSort, dSilOccupiedSort, maxSilPoint*sizeof(uint), cudaMemcpyDeviceToHost) );
	polarOccupiedSort[maxSilPoint] = 0;
	uint* p1=polarOccupiedSort, *p2=polarOccupiedPad;
	int polarPad=0, local=0, np = (*p1)&patchMask;
	while(*p2++ = *p1++)
	{
		local++;
		if( (*p1&patchMask)!=np || local==PIXELTHREADS)
		{
			for(int i=local; i<PIXELTHREADS; i++)
				*p2++ = 0;
			polarPad++;
			local = 0;
			np = (*p1&patchMask);
		}		
	}
	if(local)
	{
		for(int i=local; i<PIXELTHREADS; i++)
			*p2++ = 0;
	}
	polarPad *= PIXELTHREADS;
	//printf("%d %d\n", silSize, polarPad);

	if(polarPad*sizeof(uint)>silSize)
		printf("Allocate more GPU memeory for dSilOccupiedSortPad!\n");


	CUDA_SAFE_CALL( cudaMemset(dSilOccupiedSortPad, 0, polarPad*sizeof(uint)) );
	CUDA_SAFE_CALL( cudaMemcpy(dSilOccupiedSortPad, polarOccupiedPad, polarPad*sizeof(uint), cudaMemcpyHostToDevice) );

#if DEBUG_BUFFERS
	printf("dSilOccupiedSortPad:\n");
	dumpBuffer(dSilOccupiedSortPad, polarPad);
#endif

	Db = dim3(PIXELTHREADS, 1, 1);
	Dg = dim3(polarPad/PIXELTHREADS, 1, 1);
	//Dg = dim3(iter+1, 1, 1);
	//Dg = dim3(nPatchEveryBlock*maxPolarIndexBlock, 1, 1);

#ifdef NORMALFUNCTION
	size_t Ns = PIXELTHREADS*sizeof(float)*2;
#else
	//512&1024 work fine with deg.xyz=3;
	size_t Ns = /*(deg.x+1)*degShift.x*3*sizeof(float)*2*/1024;
#endif





	//printf("%d, %d\n", sizeof(MULTI_SIL_RENDER)/(PIXELTHREADS*sizeof(float)), sizeof(MULTI_SIL_RENDER)%(PIXELTHREADS*sizeof(float)));

#if	UNIFORM_KNOT	
	UpdatePointKernel<<<Dg, Db, Ns>>>((unsigned char*)surfaceIn, (unsigned char*)surfaceOut, pitch,
		//renderWidth, renderHeight, eye, errorXTol,
		//viewProj0, viewProj1, viewProj2, viewProj3,
		//invViewProj0, invViewProj1, invViewProj2, invViewProj3,
		dMultiSilRender, //sizeof(MULTI_SIL_RENDER)/(PIXELTHREADS*sizeof(float)), sizeof(MULTI_SIL_RENDER)%(PIXELTHREADS*sizeof(float)),
		pixelBits, patchBits, 
		//dSilOccupiedSort, dSilOccupiedSortIndex,
		dSilOccupiedSortPad, dFuncData, dPatchIndex,
		deg, degShift, seg, segShift,
		boxSize, box0, iter);
#else
#endif

	CUDA_SAFE_CALL(cudaThreadSynchronize());
	CUT_SAFE_CALL(cutStopTimer(hTimer));
	timeValue = cutGetTimerValue(hTimer);
	printf("MulitSilRender time ...: %f msec.\n", timeValue);
 

#if DEBUG_BUFFERS
	printf("dSilOccupiedSortPad:\n");
	dumpBuffer(dSilOccupiedSortPad, iter+1*64);
#endif
	//dumpBuffer(dSilOccupiedSortPad, (iter+1)*64);

}

extern "C"
void CleanUpMultiRenderCUDA()
{
	CUDA_SAFE_CALL( cudaFree(dFuncData) );
	CUDA_SAFE_CALL( cudaFree(dPatchIndex) );

	free(polarOccupiedPad);
	free(polarOccupiedSort);

	CUDA_SAFE_CALL( cudaFree(dMaxSilPoint) );
	CUDA_SAFE_CALL( cudaFree(dSilOccupied) );
	CUDA_SAFE_CALL( cudaFree(dSilOccupiedCompact) );
	CUDA_SAFE_CALL( cudaFree(dSilOccupiedValid) );
	CUDA_SAFE_CALL( cudaFree(dSilOccupiedSort) );
	CUDA_SAFE_CALL( cudaFree(dSilOccupiedSortPad) );
	//CUDA_SAFE_CALL( cudaFree(dSilOccupiedSortIndex) );

	CUDA_SAFE_CALL( cudaFree(dMultiSilRender) );

#ifndef UNIFORM_KNOT	
	CUDA_SAFE_CALL( cudaFree(dBoxX) );
	CUDA_SAFE_CALL( cudaFree(dBoxY) );
	CUDA_SAFE_CALL( cudaFree(dBoxZ) );
#endif
}
